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Abstract - We perform a matrix product state based density matrix renormalisation group 
analysis of the phases for the disordered one-dimensional Bose-Hubbard model. For particle 
densities N/L — 1, 1/2 and 2 we show that it is possible to obtain a full phase diagram using only 
the entanglement properties, which come for free when performing an update. We confirm the 
presence of Mott insulating, superfluid and Bose glass phases when N/L — 1 and 1/2 (without 
the Mott insulator) as found in previous studies. For the N/L — 2 system we find a double lobed 
superfluid phase with possible reentrance. 


Introduction. — The study of bosons in one dimen¬ 
sion has been of great interest in both theoretical and 
experimental physics for many years due in part to the 
existence of a quantum phase transition from a superfluid 
to insulator at zero temperature [l]. The introduction of 
disorder causes a further phase transition into a localised 
Bose glass phase, which is insulating but remains com¬ 
pressible [2]. The experimental study of phase transitions 
in bosonic systems is possible using Helium in porous me¬ 
dia |3 ,[4], Josephson junction arrays [ 5 ], thin films [6,7 
and, more recently, optical lattices |l ,[8]. It is now possi¬ 
ble to introduce disorder in a controlled manner to optical 
lattices using speckle potentials as to study these tran¬ 


sitions directly 11 12 


the maximum number of bosons per site is set to two. For 
disordered systems Giamarchi and Schulz used renormal¬ 
ization group (RG) techniques to determine the weak dis¬ 
order physics given the Luttinger parameter K 14,15 


some of the most effective means of garnering informa¬ 
tion. Quantum Monte Carlo has been employed in 1, 2 
and 3 dimensions 18 -23 , but these methods become diffi- 


to a number of physical systems from quantum chemistry 
to quantum information 26 , including the disordered 

The phase diagrams ob- 


25 


Bose-Hubbard model 27 28 


tained using these methods for the one dimensional case 
19 28 , whilst qualitatively agreeing, are quantitatively 
quite different. This difference could be down to the choice 
of different observables and the difficulties each has with 
finite size effects. 

In recent years the use of entanglement properties as 
a means of deciphering phase has become commonplace 
Entanglement is a measurement of a wavefunc- 


29-34 


Analytical results even for clean systems are limited. 
There is an approximate Bethe-ansatz solution 13 , where 


tion’s non-locality and as such it is an ideal means of 
analysing various phases. Modern numerical techniques 
such as tensor networks and DMRG obtain entanglement 
information as part of the update algorithms, so large 
amounts of information about the phase is gathered auto¬ 
matically in the course of the RG iterations 35 . In this 


paper we perform a DMRG simulation of the disordered 
Bose-Hubbard model in the form of a variational update 


There are further real-space RG results for the case of 
strong disorder 16 17 . Numerical approaches provide 


of a matrix product state (MPS) [35 36 implemented in 
the ITensor libraries p7j. The disordered Bose-Hubbard 
model is made up of bosonic creation, b \, and annihilation 
operators, bi , on sites of a linear lattice. The Hamiltonian 


28 


cult in the limit of zero temperature. An ideal method for 
analysing one dimensional systems is the density-matrix 
RG (DMRG) [24]. It has been applied with great success 


t ^ ^ JJ 

H = ~ X 7;( b i b i+i+h.c.)+'^2—n i (n i -l)+iJ, i n i , (1) 
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where fa = b\bi is the local occupation or number oper¬ 
ator that gives the number of bosons on site i. The po¬ 
tential disorder is modelled via uniformly distributed ran¬ 
dom chemical potentials fa G [—A/i/ 2 , A/i/ 2 ]. For ease 


of comparison, we have adopted prior conventions 28 for 


states 35 


Eg can be found by calculating the energies En+ 1 , E 


'N 


and En-i of the N + 1 , 
respectively, as 


Eg = E n +i 


2 E n + E N -\. 
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hopping t and interaction U and throughout the rest of 
the analysis t = 1 . 

Observables. — The Mott insulator can be differen¬ 
tiated from the Bose glass phase by the existence of the 
Mott gap , Eg , between the ground and first excited state. 
While DMRG ordinarily finds the ground state of the sys¬ 
tem, low lying excited states have to be constructed iter¬ 
atively by orthogonalising with respect to the lower lying 
For the Bose-Hubbard chain it is numerically 


where K is the Luttinger parameter. K takes the value 2 
for a Kosterlitz-Thouless (KT) transition from superfluid 
to Mott insulator 40 IEL By utilizing an RG approach, 
Giamarchi and Schulz 14] showed that disorder scales to 
zero in the weak disorder regime when K > 3/2, giving a 
superfluid phase. On the other hand, disorder grows for 
K < 3/2 signifying a Bose glass. This was later extended 
[15] to the medium disorder case (U ~ A/x). Instead of 
the infinite-system size result (4| we use the conformal 
field theory (CFT) expression 
size L, 


38 for an open chain of 


( b l b j) 


more convenient to use the fact that the energy of the 
excited state is equal to the difference in energy between 
the chemical potential for particle, (i v — En+i — En, and 
hole, fih = En — En- 1 , excitations 28]. This means that 


7r \/i sin (T)l i sin (y) 

2 L • • Tr(i-j) 
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1/2 K 
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N and N — 1 particle sectors, 


The expression has to be averaged over all i and j with 
separation \i—j \ and, of course, also averaged over disorder 
realisations. Calculating correlation functions does not 
require multiple DMRG runs, but requires the calculation 
of an expectation value for each combination of i and j, of 
which there are L(L — l)/2. Furthermore, the accuracy of 
locating the KT transition from correlation functions for 
the Bose-Hubbard model has previously been questioned 


Hence the determination of E g requires a DMRG run for 
each of the three different particle numbers and each set 
of parameters. 

The superfluid phase is determined by a non-zero super¬ 
fluid fraction p s . This is defined as the difference between 
the ground state energies of a chain with periodic bound¬ 
aries and anti-periodic boundaries, 


38,42 


In each DMRG run, bipartitioning the chain into sys¬ 
tem and environment blocks is done routinely to compute 
singular-value decompositions 35 . These singular values, 


s a , can themselves be used to obtain information regard¬ 


ing the phase 31 33 without the need for multiple DMRG 
runs, thus saving substantial numerical costs. The most 
common such measure is the entanglement entropy or Von 
Neumann entropy defined as 


where L is the chain length and N the number of bosons 


Sa\b = -Trp A log 2 Pa = - Yi s l lo S 2 


a= 1 


(6) 


For Mott insulator and Bose glass phases, we have 
p s = 0 , so a finite p s indicates superfluidity in the phase 
diagram. From a computational point of view, p s is not an 
easy quantity to determine as it requires the use of periodic 
boundaries, which are well-known to converge slower and 
be less accurate than for open systems when using DMRG 


which gives the entanglement between regions A and B 


Furthermore, as p s is the difference between two 
energies, two such periodic DMRG calculations have to 
be performed for each set of parameters. 

The two-point correlation function (b\bj) provides in¬ 
formation regarding the localisation of the wavefunction. 
For the Bose glass and Mott insulating phases the corre¬ 
lation function decays exponentially, ((b\bfl) oc 
where £ is the correlation length and ((...)) denotes the 
expectation value when averaged over all pairs of sites sep¬ 
arated by \i—j\ and all disorder realisations [38]. Extended 
phases like the superfluid are not localised so £ diverges in 
the thermodynamic limit. In the absence of disorder the 
superfluid phase will be described by Luttinger liquid the- 
hence the correlation function will admit a power 


35 . The reduced density matrix, pa, for region A is ob¬ 
tained from the density matrix by tracing over degrees of 
freedom from region B. Its eigenvalues are given as squares 
of the Sa s. Hence 5 a|b is a measure of the spread of the 
s a values. If there is one non-zero singular value then the 
regions are in a product state of the two regions. The 
other extreme is if all singular values are equal, in which 
case the subsystems are maximally entangled. In the sub¬ 
sequent analysis we shall average the entanglement en¬ 
tropy over all possible bipartitions along the chain. This 
averaged entanglement entropy can distinguish between 
phases with high and low entanglement, for example the 
superfluid and Mott insulating phases. 


Deng et. al. 33 used the entanglement spectral parame¬ 


ter , £, to obtain the phase diagram for an extended Bose- 
Hubbard model. The £ parameter is defined as the sum of 
the difference between the first and second, and third and 
fourth, respectively, eigenvalues s 2 a of pA when averaged 
over all bipartition positions such that La + Lb = T, i.e. 


£ — Ai — A 2 + A 3 — A 4 , 


( 7 ) 
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Fig. 1: (Color online) (Top) the largest four singular values, si (black o), S 2 (red □), s 3 (green +), and S 4 (blue x) and 
(Bottom) the entanglement entropy Sab (black A) for all possible bipartition positions, La, along a chain of length L — 50 
for (a) superfluid with U = 0.5, A fi = 1, (b) Mott insulator with U = 4.5, Ay — 1, and (c) Bose glass with U = 4.5, Ay = 7. 
The dashed horizontal line in the top (bottom) graph shows the average value of si (5a|b) while the grey shading indicates its 
standard deviation when averaged over all La positions. Solid lines connecting symbols are guides to the eye. 


with X a = Ef A i 1 Sa( I 'A)/(X- 1 ), a = 1 , 2 ,3,4, the 
bipartition-averaged a-th eigenvalue. In fig. [l] we show 
the typical behaviour of the four lowest s a values in the 
superfluid, fig. [IJa), Mott insulator, fig. [ljb) , and Bose 
glass, fig. 0 c), regimes We see that the entanglement spec¬ 
trum of the superfluid phase is somewhat noisy with the 
four singular values being of the same order of magnitude. 
Therefore the resulting 5 a |b is large, ( is small, but neither 
have too much variation along La- One of the striking fea¬ 
tures of the entanglement spectrum for the Mott insulator 
regime is that si ~ 1 while s 2 ~ S 3 ~ 54 ~ 0 for all bipar¬ 
titions, even in the presence of disorder. This means that 
the average 5 a |b will be low and ( « 1 , with a negligible 
deviation. Last, entanglement spectra of the Bose-glass 
show pronounced localised regions separated by areas of 
low entanglement. This results in a much larger variation 
of ( and 5 a|b than in the superfluid phase, with ( and 
average 5 a|b? between the other phases. These findings 
suggest that the average and spatial variations of 5a|b 
and £ might also be used to distinguish the phases of the 
disordered Bose-Hubbard model. 


Results. — We use 5a|b and C 1° create qualitative 
phase diagrams for a modest size of L = 50, disorder- 
averaged over 100 samples using as our DMRG implemen¬ 
tation the ITensor libraries 37 . The finite-size scaling 
(FSS) behaviour of 5 a|b and C I s currently not well un¬ 
derstood, thus to find the phase boundaries for L 00 , 
we perform scaling with K and E g from estimates up to 
L = 200. This is a numerically more expensive proce¬ 
dure, so we concentrate on a small number of points with 
positions motivated by the phase diagram from the entan¬ 
glement properties. For disordered systems, getting stuck 
in local minima is particularly problematic, so we use a 
relatively large bond dimension x = 200 and perform 20 
DMRG sweeps of the chain for each sample. Our trunca¬ 
tion error is less than 10 -10 . We also introduce a small 
noise term for the first few sweeps; this perturbs a perhaps 
bad initial wavefunction, allowing faster convergence into 
the ground state. 

Bosons do not obey the Pauli exclusion principle and 


hence can condense onto a single site. In order to capture 
such a behaviour, the one-site basis dimension should to 
be as large as the number of particles in the system. This 
is numerically infeasible and it is necessary to introduce a 
finite maximum number of bosons that can occupy each 
site. We use max(n^) = 5, consistent with ref. 41 who 


find that a higher particle number does not effect the re¬ 
sults appreciably for U > 0 SUM 

Density = 1. For particle density N/L = 1 , in the 
clean case, the system is in a superfluid phase for small 
U but transitions into a Mott insulating phase at a crit¬ 
ical U c • Introducing disorder enables the existence of a 
localized Bose glass phase [2 . The possibility of a direct 
transition from superfluid to Mott insulator has been dis¬ 
cussed extensively (see references in [43]). In one dimen¬ 
sion it was shown [44] that the transition necessarily goes 
via the Bose glass phase. This is now also the accepted 
picture for any dimension 


43 


We show our results based on £ and 5 a |b f° r L = 50 
in fig. [ 2 ] (a) and (b), respectively. The superfluid, small 
U <1.5, and the Mott insulator, U > 2, are clearly distin¬ 
guishable in both panels. The boundary of the superfluid 
to the Bose glass is less well defined and it is not clear 
that there is a Bose glass region between the Mott insu¬ 
lator and the superfluid — very different wavefunctions 
give similar average entanglement entropy. Following on 
from our prior discussion of fig. [l] we also plot in fig. [2jc) 
the standard erro]0 of £, A£, and, similarly, (d) A5a|b- 
In these plots the phases become clear and their bound¬ 
aries are consistent with earlier work 28 . In particular, 


a Bose glass phase can be easily identified between Mott 
insulator and superfluid. Furthermore, we see that the 
contours for £ and 5a|b in fig- [2] are qualitatively simi¬ 
lar, just as those for AC and A5 a| B . We emphasize that 
for the entanglement-based measures present here, it is in 
fact possible to discern all of the phases with just a sin- 


x Five is the current hard limit in the ITensor code. 

2 We note that for larger system sizes the variance or standard 
deviation may be better measures of distribution width as they do 
not approach zero in the infinite system limit. 
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Fig. 2: (Color online) Phase diagrams for the disordered Bose-Hubbard model at N/L — 1 given as contour plots of (a) £, (b) 
Sa|b> ( c ) A£ and (d) ASa\b- The color shading goes from low (orange/dark) to high (blue/white) value and its coarse-graining 
reflects the (U, A/i) resolution of our calculations for L — 50. The contour lines correspond to (a) £ =s 0.1, 0.2,..., 0.9, (b) 
Sa\b — 1.4,..., 0.4,0.3, (c) A£/10 -4 = 0 .1, 1 , 2 ,..., 10 and (d) A*Sa|b/10 -4 = 0 .1, 1 , 2 ,..., 10 . In all cases the two extreme 
contours values are shown as dashed lines. Note that regions of high £ corresponds to low Sa\b- The circles (white) and squares 
(blue) denote estimations of K and E g from FSS for L —>• oo while the stars (red) indicate the Km 2 values for L = 50 as 
discussed in the text. The arrow (black) denotes the expected transition in the clean case at U c - The dotted straight line 
indicates Ap = 2 U. Error bars (white) show the standard error of the mean in all cases. They are within symbol size if not 
shown. We emphasize that the color shading does not directly indicate the transitions, but rather quantifies the change in 
entanglement measures. The grey line highlights the start of the shaded region where the probability of (rii) > 4.9 is greater 
than 10 -3 . 


gle DMRG run for each (U, A/i, disorder realisation) data 
point. This is a clear advantage in terms of numerical 
costs when compared to calculations based on E gi p s or 
K. 

In order to augment the finite-size phases identified in 
fig- HI we now perform runs with larger L and employ 
FSS. To find the superfluid-Bose glass transition in the 
thermodynamic limit we calculate K for various points 
along the boundary for system sizes L = 30, 50, 100 and 
150. The transition is of KT type at K = 3/2. The 
corresponding points in ([/, A/i) which are shown as filled 
circles in fig. [2j For reference, we also plot the points where 
K = 3/2 for L = 50 (stars). Similarly, the superfluid- 
Mott insulator transition point U c is the point on the zero 
disorder axis where K = 2. We estimate it value as U c = 
1.634 =b 0.002. We also calculate E g for the same system 
sizes and use FSS to find the Mott insulator-Bose glass 
boundary indicated as squares in fig. |2j The superfluid 
region we find is significantly smaller than that of ref. 28 
but matches ref. 19 ; the position of Mott insulator-Bose 


glass boundary is very similar to 28 and different to 
The RG analysis of Refs. 


19 


14 


and 15 suggests that 


there may be a further Anderson glass phase in the low 
U < A/i/2 region of the Bose glass phase highlighted by 
the dashed line in fig. [2j This would imply a critical point 
along the superfluid boundary at which point K at the 
transition becomes disorder dependent. Our entanglement 
analysis shows no sign of such a transition either within 
the Bose glass phase or on the boundary with the super¬ 
fluid. However, when U <C A/i the truncation of the basis, 
i.e. max(n^) < 5, becomes more problematic so we cannot 
rule out the existence of another phase in this region. 

Density = 1/2. The clean case for N/L = 1/2 remains 
a superfluid for all values of U [2]. When A/i is increased, 


our entanglement measures indicate the eventual emer¬ 
gence of a Bose glass phase as shown in fig. |3|a+b). Still, 
the superfluid phase for L = 50 seems to extend up to 
A/i < 1 for U < 5 as shown by all four entanglement mea¬ 
sures. The Giamarchi-Schulz criterion 14p8 implies that 
the Bose-Hubbard model should be in a Bose glass phase 
for K < 3/2. In fig. [^a-fb) we show that the resulting 
boundaries indicate that the superfluid phase extends as 
far as Uk= 3/2 = 3.5 ±0.1, i.e. it ends somewhat earlier 
for low A/i than suggested by our entanglement measures. 
In order to explore this region further, we have also cal¬ 
culated p s for fixed A/i = 0.5 and sizes L = 50, 100, 150, 
and 200 as shown in fig. Qa) . The results for p s have been 
computed for increased bond dimension x = 400 with 40 
DMRG sweeps and 20 disorder configuration to offset the 
reduction in precision due to periodic boundaries. The 
figure shows that for U > 3, p s decreases when increas¬ 
ing L as expected in the Bose glass phase. However, the 
decrease is very slow and, for the system sizes attainable 
by us, even seems to saturate at non-zero values. These 
results suggest that for finite systems, the K = 3/2 crite¬ 
rion significantly underestimates the extent of the super¬ 
fluid phase, while our four entanglement measures and p s 
predict a much larger region. Performing a FSS analysis 
for K as shown in fig. [4^a) we find the U values, for which 
K = 3/2 in the limit L —>• 00, converge towards a limit¬ 
ing value of U c = 3.09 ± 0.01 (see also fig. [3|. This again 
indicates that in an infinite system, we expect the super¬ 
fluid to Bose glass transition to take place at much lower 
values of U c than observed for L = 50. The relevance of 
this result is of course that experimental realizations of the 
Bose-Hubbard model are typically in cold atom systems, 
which are limited to finite system sizes, currently a typical 
lattice dimension is ^ 50 — 100 [l]. 
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Fig. 3: (Color online) Phase diagrams for the disordered Bose-Hubbard model at (a+b) N/L — 1/2 and (c+d) N/L — 2 given 
as contour plots of (a+c) £ and Sa\b, (b+d) A£ and ASa\b • Colors, symbols and lines (solid and dashed) denote corresponding 
estimates as in fig. i The black contour lines correspond to (a) £ = 0.2, 0.3,..., 0.8, (b) A£/10 3 = 0.1,1,1.5, 2, 3,4 for 
N/L = 1/2 and (c) £ = 0.1, 0.2,..., 0.9, (d) A£/10 -4 = 0.1,1, 2, 3,4, 6, 8, 9,10 for N/L = 2. The white contour lines represent 
the results for (a+c) Sa\b, (b+d) ASa\b- For Sa\b , the contour values are 1.2,..., 0.4, 0.3 in (a), while they are 1.4,..., 0.3, 0.2 
for (c); for ASa\b the values are as for A£. High £ corresponds to low Sa\b- The black arrow corresponds to U c for both 
densities as discussed in the text. The two dotted straight lines indicates A p — 2 U and 4 U. The grey line and area have the 
same meaning as in fig. [2] 


For values of A/i > 1, the situation is less severe and we 
see in fig. [3ja+b) that our entanglement-based measures 
again qualitatively agree with the Giamarchi-Schulz crite¬ 
rion, both for L = 50 and estimated via FSS at L -+ oo. 

Density = 2. To the best of our knowledge, the phases 
for N/L = 2 have not been shown before in the literature. 
Due to our numerical restriction of five bosons per site, 
this regime is close to the limit of what can be studied 
reliably, particularly for small U where the occupancy per 
site should be large. For large J7, one might expect that 
we will have a Mott insulator of boson pairs , while a super¬ 
fluid of boson pairs emerges for small U and small A/i. As 
before, we envisage a disordered Bose glass phase for large 
A/i. With more particles per site than in the N/L = 1 
case, we could furthermore expect that onset of the Mott 
transition at A/i = 0 is at larger values of £/, since there 
is a larger energy penalty to pay for a doubly occupied 
site |2 . Similarly, as the cost for two boson pairs to go 
onto the same site is 2U, we expect the 2 U = A/i/2 line to 
characterize the superfluid phase as in the N/L = 1 case. 
In addition, one might conjecture to see a remnant of the 
U = A/i/2 condition. 

In fig.^c+d), we show that our expectations are largely 
validated. In particular, a double lobe shape for the su¬ 
perfluid phase emerges and allows a possible re-entrant 
behaviour given a suitable cut across parameter space. 
The gradient of the Mott insulating phase boundary is 
shallower (^ 4/3) when compared to N/L = 1. Further¬ 
more, both the £ and 5 a|b based entanglement measures, 
as well as their errors, A£ and A5 a|b> capture the phases 
equally well and agree with the K and E g estimates. Note 
that for N/L = 2, the KT superfluid-to-Mott transition 
at A/i = 0 corresponds to K = 2 and we finite-size scale 
the Luttinger parameter to find U c = 2.75 =b 0.03. 

We emphasize that the points for small I/, see top left 
of the phase diagram in fig.[3jc+d), should be viewed with 



Fig. 4: (Color online) (a) Superfluid fraction p s (JJ ) for N/L — 
1/2 with A/i = 0.5 for lengths 50-200. The vertical line 
indicates U c = 3.09. The inverted triangles give the finite- 
size scaled L oo limit with the dashed line a guide to the 
eye. (b) Luttinger parameter K for various lengths 30-150 at 
N/L — 1/2. The horizontal line highlights K — 3/2. The inset 
shows the FSS analysis. 


caution as the basis truncation will affect the results. The 
grey line in fig. [3] — as in fig. [2]— indicates the points at 
which the probability of obtaining a site with (rii) > 4.9 
reaches 10 -3 . This clearly shows that for the Bose glass 
with small V all wavefunctions are beginning to reach the 
limit five of bosons per site, however in the bulk of the 
phase diagram the results are not affected. 


Conclusion. — We have analysed the phase diagrams 
of the disordered Bose-Hubbard model for fillings N/L = 
1/2, 1 and 2 using the entanglement-based measures £, 
5 a|b> A£ and A5a|b- We find that despite success in 
ref. 


33 


£ or 5 a|b alone do not always faithfully repro¬ 
duce" tSe phase diagrams. The error-based measures, A£ 
and A5 a|b? provide a much clearer picture — the distri¬ 
butions of the values contain more information regarding 
the nature of the phase than the mean values alone. These 
measures are an excellent means of quickly identifying the 
different phases of the system while removing the need 
for multiple DMRG runs per measurement and special 
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boundary conditions. Unfortunately, they do not seem to 
exhibit a simple FSS behavior, at least for the system size 
up to L = 200 used here. While ( and 5a|b and, in par¬ 
ticular, and A5 a|b provide a numerically convenient, 
qualitative outline of the phase boundaries, it seems still 
necessary to apply FSS to K and E g for estimates of the 
boundaries in the L -A oo limit. For N/L = 1 our phase 
diagram is found to complement the results of refs. 19,28 


For N/L = 1/2 the diagram shows strong finite-size effects 
and the critical U defined by the Giamarchi-Schulz crite¬ 
rion is not apparent for these finite systems. Finally, for 
N/L = 2 the superfluid phase has a double-lobed appear¬ 
ance giving rise to re-entrance phenomena. 
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